IODS course project
Jaakko Keinanen
About the project
This seems like a great course, although I don??t know yet whether I will have immediate use for the skills I learn here but hopefully in the future this will be useful. I??m also sure that learning about R and open data will make me a better researcher.
http://www.github.com/jaakkoik/IODS-project
1. Description of the data
The dataset consists of students responses to a questionnaire about study habits and attitudes, including also exam points reflecting how succesful they were during a course. It consists of 8 variables, including the id-variable X, gender, age, and variables that describe the students attitude and learning styles, and exam points. The data has 166 different students answers in it.
list.files("~/Documents/Jaakon tutkimus/Open data kurssi s2017/GitHub/IODS-project/data")## [1] "create_alc.R" "create_alc.csv" "learning2014.csv"
## [4] "student-mat.csv" "student-por.csv"students2014 <- read.csv("~/Documents/Jaakon tutkimus/Open data kurssi s2017/GitHub/IODS-project/data/learning2014.csv", header =TRUE, sep = ",")
dim(students2014)## [1] 166 8str(students2014)## 'data.frame': 166 obs. of 8 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ gender : Factor w/ 2 levels "F","M": 1 2 1 2 2 1 2 1 2 1 ...
## $ Age : int 53 55 49 53 49 38 50 37 37 42 ...
## $ Attitude: int 37 31 25 35 37 38 35 29 38 21 ...
## $ deep : num 3.58 2.92 3.5 3.5 3.67 ...
## $ stra : num 3.38 2.75 3.62 3.12 3.62 ...
## $ surf : num 2.58 3.17 2.25 2.25 2.83 ...
## $ Points : int 25 12 24 10 22 21 21 31 24 26 ...2. Graphical overview of the data:
library(ggplot2)
library(GGally)
p3 <- ggpairs(students2014, mapping = aes(col=gender, alpha=0.3), lower = list(combo = wrap("facethist", bins = 20)))
p3
From the graph we see that exam points correlates most strongly with attitude and the second largest correlation is with strategic learning -variable. From the distributions it is visible that there are some differences between sexes in the distributions, for example men tend to have more attitude points but the differences in most of the variables are not large.
3. The regression model:
Dependent variable: points. Independent variable: attitude, gender, stra.
my_model <- lm(Points ~ Attitude + gender + stra, data = students2014)
summary(my_model)##
## Call:
## lm(formula = Points ~ Attitude + gender + stra, data = students2014)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.7179 -3.3285 0.5343 3.7412 10.9007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.97982 2.40303 3.737 0.000258 ***
## Attitude 0.35100 0.05956 5.893 2.13e-08 ***
## genderM -0.22362 0.92477 -0.242 0.809231
## stra 0.89107 0.54409 1.638 0.103419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.304 on 162 degrees of freedom
## Multiple R-squared: 0.2051, Adjusted R-squared: 0.1904
## F-statistic: 13.93 on 3 and 162 DF, p-value: 3.982e-08From the summary of the model we see that only Attitude is significant (Pr < 0.05). Gender is the least significant (Pr largest). We run the model again after removing gender from it.
my_model2 <- lm(Points ~ Attitude + stra, data = students2014)
summary(my_model2) ##
## Call:
## lm(formula = Points ~ Attitude + stra, data = students2014)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.6436 -3.3113 0.5575 3.7928 10.9295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.97290 2.39591 3.745 0.00025 ***
## Attitude 0.34658 0.05652 6.132 6.31e-09 ***
## stra 0.91365 0.53447 1.709 0.08927 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.289 on 163 degrees of freedom
## Multiple R-squared: 0.2048, Adjusted R-squared: 0.1951
## F-statistic: 20.99 on 2 and 163 DF, p-value: 7.734e-09After removing gender from the model stra is still not significant but the coefficient estimate is slightly larger.
- From the summary of the last model (my_model2) we can interpret that when Attitude points increase by one the exam points increase by 0.35 (the coeffficient estimate is 0.35). When stra points increase by one the exam points increase by 0.91. The multiple R squared is 0.20 which shows that 20 % of the variability in the exam scores is explained by the attitude and stra variables.
5. Diagnostic plots.
plot(my_model2, which=c(1, 2, 5))
QQ-plot shows how the errors in the model are distributed. The errors should be normally distributed for the regression model to be reliable. From the plot it is visible that the circles fit on the line quite well although there is some deviation from the line in both ends. All in all, the errors seem to be normally distributed.
Residual vs fitted -plot shows how well the data meets the assumption that the size of the errors should not depend on the explanatory variables in the model. When the points are randomly spread on the plot, as they seem to be here, the assumption is usually met
Residuals vs leverage plot shows how big of an impact the outliers get. In the plot, circles in the upper right or lower left corners may be influential outliers. When circles are out of cooks distance the are influential to the regression model. Here all the cases are within cook??s distance so the assumption is met.
Chapter 3. RStudio Exercise 3
Read data to “alc” and check the column names.
alc <- read.csv("/Users/caroram/Documents/Jaakon tutkimus/Open data kurssi s2017/GitHub/IODS-project/data/create_alc.csv")
colnames(alc)## [1] "school" "sex" "age" "address" "famsize"
## [6] "Pstatus" "Medu" "Fedu" "Mjob" "Fjob"
## [11] "reason" "nursery" "internet" "guardian" "traveltime"
## [16] "studytime" "failures" "schoolsup" "famsup" "paid"
## [21] "activities" "higher" "romantic" "famrel" "freetime"
## [26] "goout" "Dalc" "Walc" "health" "absences"
## [31] "G1" "G2" "G3" "alc_use" "high_use"str(alc)## 'data.frame': 382 obs. of 35 variables:
## $ school : Factor w/ 2 levels "GP","MS": 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "F","M": 1 1 1 1 1 2 2 1 2 2 ...
## $ age : int 18 17 15 15 16 16 16 17 15 15 ...
## $ address : Factor w/ 2 levels "R","U": 2 2 2 2 2 2 2 2 2 2 ...
## $ famsize : Factor w/ 2 levels "GT3","LE3": 1 1 2 1 1 2 2 1 2 1 ...
## $ Pstatus : Factor w/ 2 levels "A","T": 1 2 2 2 2 2 2 1 1 2 ...
## $ Medu : int 4 1 1 4 3 4 2 4 3 3 ...
## $ Fedu : int 4 1 1 2 3 3 2 4 2 4 ...
## $ Mjob : Factor w/ 5 levels "at_home","health",..: 1 1 1 2 3 4 3 3 4 3 ...
## $ Fjob : Factor w/ 5 levels "at_home","health",..: 5 3 3 4 3 3 3 5 3 3 ...
## $ reason : Factor w/ 4 levels "course","home",..: 1 1 3 2 2 4 2 2 2 2 ...
## $ nursery : Factor w/ 2 levels "no","yes": 2 1 2 2 2 2 2 2 2 2 ...
## $ internet : Factor w/ 2 levels "no","yes": 1 2 2 2 1 2 2 1 2 2 ...
## $ guardian : Factor w/ 3 levels "father","mother",..: 2 1 2 2 1 2 2 2 2 2 ...
## $ traveltime: int 2 1 1 1 1 1 1 2 1 1 ...
## $ studytime : int 2 2 2 3 2 2 2 2 2 2 ...
## $ failures : int 0 0 2 0 0 0 0 0 0 0 ...
## $ schoolsup : Factor w/ 2 levels "no","yes": 2 1 2 1 1 1 1 2 1 1 ...
## $ famsup : Factor w/ 2 levels "no","yes": 1 2 1 2 2 2 1 2 2 2 ...
## $ paid : Factor w/ 2 levels "no","yes": 1 1 2 2 2 2 1 1 2 2 ...
## $ activities: Factor w/ 2 levels "no","yes": 1 1 1 2 1 2 1 1 1 2 ...
## $ higher : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ romantic : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 1 1 1 ...
## $ famrel : int 4 5 4 3 4 5 4 4 4 5 ...
## $ freetime : int 3 3 3 2 3 4 4 1 2 5 ...
## $ goout : int 4 3 2 2 2 2 4 4 2 1 ...
## $ Dalc : int 1 1 2 1 1 1 1 1 1 1 ...
## $ Walc : int 1 1 3 1 2 2 1 1 1 1 ...
## $ health : int 3 3 3 5 5 5 3 1 1 5 ...
## $ absences : int 5 3 8 1 2 8 0 4 0 0 ...
## $ G1 : int 2 7 10 14 8 14 12 8 16 13 ...
## $ G2 : int 8 8 10 14 12 14 12 9 17 14 ...
## $ G3 : int 8 8 11 14 12 14 12 10 18 14 ...
## $ alc_use : num 1 1 2.5 1 1.5 1.5 1 1 1 1 ...
## $ high_use : logi FALSE FALSE TRUE FALSE FALSE FALSE ...Four variables that could be associated with alcohol use and hypotheses: 1) sex: boys might drink more alcohol tha girls 2) age: those who are old enough might drink more 3) studytime: those who study more might dring less 4) famrel: those with good family relationships might drink more Check distributions and correlations.
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:GGally':
##
## nasa## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(tidyr)
g1 <- ggplot(data = alc, aes(x = alc_use)) + geom_bar(aes(fill=sex), position="dodge")
g1
g2 <- ggplot(data = alc, aes(x = studytime)) + geom_bar()
g2
g3 <- ggplot(alc, aes(studytime, alc_use)) + geom_boxplot(aes(group=studytime))
g3
#logistic regression to predict alcohol high use with preselected 4 variables:
First change the studytime and famrel variables as factors (categorical).
alc$studytime <- factor(alc$studytime)
alc$famrel <- factor(alc$famrel)
m <- glm(high_use ~ age + studytime + sex + famrel, data = alc, family = "binomial")
summary(m)##
## Call:
## glm(formula = high_use ~ age + studytime + sex + famrel, family = "binomial",
## data = alc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4630 -0.8367 -0.6373 1.1158 2.2807
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.92781 1.90605 -2.585 0.00973 **
## age 0.23811 0.10162 2.343 0.01912 *
## studytime2 -0.26540 0.27519 -0.964 0.33484
## studytime3 -1.18329 0.45610 -2.594 0.00948 **
## studytime4 -1.11392 0.59560 -1.870 0.06145 .
## sexM 0.73609 0.25197 2.921 0.00348 **
## famrel2 0.75910 0.97477 0.779 0.43613
## famrel3 0.55588 0.89199 0.623 0.53316
## famrel4 0.01567 0.86874 0.018 0.98561
## famrel5 -0.29474 0.88756 -0.332 0.73983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 465.68 on 381 degrees of freedom
## Residual deviance: 426.33 on 372 degrees of freedom
## AIC: 446.33
##
## Number of Fisher Scoring iterations: 4Interpretation: Age is a significant predictor of high use (p<0.05). So is studytime category 3 (compared to category 1 which is not shown). This meand studying 5-10 hours weekly is associated with less high use than studying <2 hours. Category 4 (studying more than 10 hours weekly) was not significantly associated with high use. Male sex increased the probability of high use significantly. Family relations are not associated with high use. #Odds ratios and confidence intervals:
OR <- coef(m) %>% exp
CI<- confint(m) %>% exp## Waiting for profiling to be done...cbind(OR, CI)## OR 2.5 % 97.5 %
## (Intercept) 0.007242324 0.0001543051 0.2828385
## age 1.268854575 1.0416646124 1.5534993
## studytime2 0.766901440 0.4478608817 1.3200182
## studytime3 0.306269593 0.1186265549 0.7220426
## studytime4 0.328270171 0.0888754833 0.9696930
## sexM 2.087754319 1.2785290293 3.4399886
## famrel2 2.136357015 0.3477827163 18.2954416
## famrel3 1.743477847 0.3408184256 13.2361471
## famrel4 1.015793214 0.2089632663 7.4639818
## famrel5 0.744723795 0.1466521034 5.6113052CI:s of predictors do not include 0.
# predict() the probability of high_use
probabilities <- predict(m, type = "response")
# add the predicted probabilities to 'alc'
alc <- mutate(alc, probability = probabilities)
# use the probabilities to make a prediction of high_use
alc <- mutate(alc, prediction = probability>0.5)
# see the last ten original classes, predicted probabilities, and class predictions
select(alc, studytime, famrel, sex, high_use, probability, prediction) %>% tail(10)## studytime famrel sex high_use probability prediction
## 373 1 4 M FALSE 0.5861623 TRUE
## 374 1 5 M TRUE 0.4500661 FALSE
## 375 3 5 F FALSE 0.1071889 FALSE
## 376 1 4 F FALSE 0.3483997 FALSE
## 377 3 5 F FALSE 0.1321975 FALSE
## 378 2 4 F FALSE 0.2908050 FALSE
## 379 2 2 F FALSE 0.4630561 FALSE
## 380 2 1 F FALSE 0.2875839 FALSE
## 381 1 2 M TRUE 0.6491547 TRUE
## 382 1 4 M TRUE 0.5274744 TRUE# tabulate the target variable versus the predictions
table(high_use = alc$high_use, prediction = alc$prediction)## prediction
## high_use FALSE TRUE
## FALSE 251 17
## TRUE 89 25# initialize a plot of 'high_use' versus 'probability' in 'alc'
g <- ggplot(alc, aes(x = probability, y = high_use, col=prediction))
# define the geom as points and draw the plot
g + geom_point()
# tabulate the target variable versus the predictions
table(high_use = alc$high_use, prediction = alc$prediction) %>% prop.table %>% addmargins()## prediction
## high_use FALSE TRUE Sum
## FALSE 0.65706806 0.04450262 0.70157068
## TRUE 0.23298429 0.06544503 0.29842932
## Sum 0.89005236 0.10994764 1.00000000# define a loss function (mean prediction error)
loss_func <- function(class, prob) {
n_wrong <- abs(class - prob) > 0.5
mean(n_wrong)
}
# call loss_func to compute the average number of wrong predictions in the (training) data
loss_func(class = alc$high_use, prob = alc$probability)## [1] 0.2774869Exercise 4 Clustering and classification
The dataset Boston has 506 observations and 14 variables. It is about housing values in sub-urbs of Boston and contains information about, for example, crime rate, nitrogen oxide concentrations, pupil teacher ratio in different towns.
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(corrplot)## corrplot 0.84 loadedlibrary(tidyverse)## -- Attaching packages -------------------------------------------- tidyverse 1.2.1 --## <U+221A> tibble 1.3.4 <U+221A> purrr 0.2.4
## <U+221A> readr 1.1.1 <U+221A> stringr 1.2.0
## <U+221A> tibble 1.3.4 <U+221A> forcats 0.2.0## -- Conflicts ----------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x MASS::select() masks dplyr::select()library(dplyr)
data("Boston")
str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...dim(Boston)## [1] 506 14pairs(Boston)
summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00#correlation matrix of variables rounded in two digits
cor_matrix<-cor(Boston) %>% round(digits=2)
# print the correlation matrix
cor_matrix## crim zn indus chas nox rm age dis rad tax
## crim 1.00 -0.20 0.41 -0.06 0.42 -0.22 0.35 -0.38 0.63 0.58
## zn -0.20 1.00 -0.53 -0.04 -0.52 0.31 -0.57 0.66 -0.31 -0.31
## indus 0.41 -0.53 1.00 0.06 0.76 -0.39 0.64 -0.71 0.60 0.72
## chas -0.06 -0.04 0.06 1.00 0.09 0.09 0.09 -0.10 -0.01 -0.04
## nox 0.42 -0.52 0.76 0.09 1.00 -0.30 0.73 -0.77 0.61 0.67
## rm -0.22 0.31 -0.39 0.09 -0.30 1.00 -0.24 0.21 -0.21 -0.29
## age 0.35 -0.57 0.64 0.09 0.73 -0.24 1.00 -0.75 0.46 0.51
## dis -0.38 0.66 -0.71 -0.10 -0.77 0.21 -0.75 1.00 -0.49 -0.53
## rad 0.63 -0.31 0.60 -0.01 0.61 -0.21 0.46 -0.49 1.00 0.91
## tax 0.58 -0.31 0.72 -0.04 0.67 -0.29 0.51 -0.53 0.91 1.00
## ptratio 0.29 -0.39 0.38 -0.12 0.19 -0.36 0.26 -0.23 0.46 0.46
## black -0.39 0.18 -0.36 0.05 -0.38 0.13 -0.27 0.29 -0.44 -0.44
## lstat 0.46 -0.41 0.60 -0.05 0.59 -0.61 0.60 -0.50 0.49 0.54
## medv -0.39 0.36 -0.48 0.18 -0.43 0.70 -0.38 0.25 -0.38 -0.47
## ptratio black lstat medv
## crim 0.29 -0.39 0.46 -0.39
## zn -0.39 0.18 -0.41 0.36
## indus 0.38 -0.36 0.60 -0.48
## chas -0.12 0.05 -0.05 0.18
## nox 0.19 -0.38 0.59 -0.43
## rm -0.36 0.13 -0.61 0.70
## age 0.26 -0.27 0.60 -0.38
## dis -0.23 0.29 -0.50 0.25
## rad 0.46 -0.44 0.49 -0.38
## tax 0.46 -0.44 0.54 -0.47
## ptratio 1.00 -0.18 0.37 -0.51
## black -0.18 1.00 -0.37 0.33
## lstat 0.37 -0.37 1.00 -0.74
## medv -0.51 0.33 -0.74 1.00# visualize the correlation matrix
corrplot(cor_matrix, method="circle", type = "upper", cl.pos="b", tl.pos="d", tl.cex=0.6)
From the above correlation matrix: home value (medv) has a strong negative correlation with lower status of the population (lstat, red circle) and a strong positive correlation with average number of rooms per dwelling (rm, blue circle).
hist(Boston$medv)
From the histogram: the distribution of housing value is right-skewed.
The scaled variables all have mean of 0.0.
boston_scaled <- scale(Boston)
summary(boston_scaled)## crim zn indus
## Min. :-0.419367 Min. :-0.48724 Min. :-1.5563
## 1st Qu.:-0.410563 1st Qu.:-0.48724 1st Qu.:-0.8668
## Median :-0.390280 Median :-0.48724 Median :-0.2109
## Mean : 0.000000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.007389 3rd Qu.: 0.04872 3rd Qu.: 1.0150
## Max. : 9.924110 Max. : 3.80047 Max. : 2.4202
## chas nox rm age
## Min. :-0.2723 Min. :-1.4644 Min. :-3.8764 Min. :-2.3331
## 1st Qu.:-0.2723 1st Qu.:-0.9121 1st Qu.:-0.5681 1st Qu.:-0.8366
## Median :-0.2723 Median :-0.1441 Median :-0.1084 Median : 0.3171
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.:-0.2723 3rd Qu.: 0.5981 3rd Qu.: 0.4823 3rd Qu.: 0.9059
## Max. : 3.6648 Max. : 2.7296 Max. : 3.5515 Max. : 1.1164
## dis rad tax ptratio
## Min. :-1.2658 Min. :-0.9819 Min. :-1.3127 Min. :-2.7047
## 1st Qu.:-0.8049 1st Qu.:-0.6373 1st Qu.:-0.7668 1st Qu.:-0.4876
## Median :-0.2790 Median :-0.5225 Median :-0.4642 Median : 0.2746
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.6617 3rd Qu.: 1.6596 3rd Qu.: 1.5294 3rd Qu.: 0.8058
## Max. : 3.9566 Max. : 1.6596 Max. : 1.7964 Max. : 1.6372
## black lstat medv
## Min. :-3.9033 Min. :-1.5296 Min. :-1.9063
## 1st Qu.: 0.2049 1st Qu.:-0.7986 1st Qu.:-0.5989
## Median : 0.3808 Median :-0.1811 Median :-0.1449
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.4332 3rd Qu.: 0.6024 3rd Qu.: 0.2683
## Max. : 0.4406 Max. : 3.5453 Max. : 2.9865# change the object to data frame
boston_scaled <- as.data.frame(boston_scaled)
# create a quantile vector of crime rate
bins <- quantile(boston_scaled$crim)
bins## 0% 25% 50% 75% 100%
## -0.419366929 -0.410563278 -0.390280295 0.007389247 9.924109610# create a categorical variable 'crime' using the quantiles as breaks
crime <- cut(boston_scaled$crim, breaks = c(bins), include.lowest = TRUE, label=c("low", "med_low", "med_high", "high"))
# remove original crim from the dataset
boston_scaled <- dplyr::select(boston_scaled, -crim)
# add the new categorical value to scaled data
boston_scaled <- data.frame(boston_scaled, crime)
# number of rows in the Boston dataset
n <- nrow(boston_scaled)
# choose randomly 80% of the rows
ind <- sample(n, size = n * 0.8)
# create train set
train <- boston_scaled[ind,]
# create test set
test <- boston_scaled[-ind,]# linear discriminant analysis
lda.fit <- lda(crime ~ ., data = train)
# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "red", tex = 0.75, choices = c(1,2)){
heads <- coef(x)
arrows(x0 = 0, y0 = 0,
x1 = myscale * heads[,choices[1]],
y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
text(myscale * heads[,choices], labels = row.names(heads),
cex = tex, col=color, pos=3)
}
# target classes as numeric
classes <- as.numeric(train$crime)
# plot the lda results
plot(lda.fit, dimen = 2, col = classes, pch=classes)
lda.arrows(lda.fit, myscale = 5)
#save crime categories from the test set
crime_cat_test <- test$crime
#remove the categorical crime variable from the test dataset
test <- dplyr::select(test, -crime)
# predict classes with test data
lda.pred <- predict(lda.fit, newdata = test)
# cross tabulate the results
t<- table(correct = crime_cat_test, predicted = lda.pred$class)
100*prop.table(t, margin=1)## predicted
## correct low med_low med_high high
## low 72.000000 28.000000 0.000000 0.000000
## med_low 15.384615 50.000000 34.615385 0.000000
## med_high 4.545455 27.272727 63.636364 4.545455
## high 0.000000 0.000000 3.448276 96.551724The model predicts correctly 52% of low cases, 63% of med_low, 58% of med_high and 100% of high cases. In every category it is better than just guessing and in high category it predicts all cases correctly.
#Reload the Boston dataset and standardize the dataset
data("Boston")
summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00boston_scaled <- scale(Boston)
#Calculate the distances between the observations
dist_eu <- dist(Boston)
summary(dist_eu)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.119 85.624 170.539 226.315 371.950 626.047# k-means clustering with 4 centers
km <-kmeans(Boston, centers = 4)
# plot the Boston dataset with clusters
pairs(Boston, col = km$cluster)
#Reducing the number of pairs in the plot
pairs(Boston[1:5], col=km$cluster)
pairs(Boston[6:10], col = km$cluster)
NA
km <-kmeans(Boston, centers = 2)
pairs(Boston[1:5], col = km$cluster)
pairs(Boston[6:10], col = km$cluster)